#################################################################################### 
####################################################################################
####################################################################################
# Main Text: Figure 3--5
####################################################################################
####################################################################################
library(here)
library(kableExtra)
library(lmtest)
library(multiwayvcov)
library(ggplot2)
library(stringr)
####################################################################################
####################################################################################
## Read in survey data (from: 8_2020_plots.R)
dat = readRDS(here("data/dat_unrestricted.rds"))
datb = readRDS(here("data/dat_restricted.rds"))

####################################################################################
# President Power Plot
d2 <- dat[!is.na(dat$power_pres),] %>% 
  group_by(treatment,power_pres) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(dat[!is.na(dat$power_pres) & dat$treatment == 0,]$power_pres)) %>%
  mutate(meant=mean(dat[!is.na(dat$power_pres) & dat$treatment == 1,]$power_pres)) %>%
  mutate(sample="Full (n=1477)")


d2b <- datb[!is.na(datb$power_pres),] %>% 
  group_by(treatment,power_pres) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(datb[!is.na(datb$power_pres) & datb$treatment == 0,]$power_pres)) %>%
  mutate(meant=mean(datb[!is.na(datb$power_pres) & datb$treatment == 1,]$power_pres)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

pres_power = ggplot(d2, aes(x = factor(power_pres), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  theme_bw() + 
  geom_text(aes(x=1.5, y=50, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=1.5, y=46, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("None","A little","Some","A lot")) + 
  labs(#title = "Perception of President's Power over Where NGOs Locate Projects",
       x = "Perceived Power of President over NGO Projects",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "black"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.1, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/figure_3.tiff"), pres_power, width = 8, height = 5, dpi = 700)
ggsave(file = here("plots/pres_power_wrap.pdf"), pres_power, width = 8, height = 5)

####################################################################################
# President health performance plot
d2 <- dat[!is.na(dat$performance_health_pres),] %>% 
  group_by(treatment,performance_health_pres) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(dat[!is.na(dat$performance_health_pres) & dat$treatment == 0,]$performance_health_pres)) %>%
  mutate(meant=mean(dat[!is.na(dat$performance_health_pres) & dat$treatment == 1,]$performance_health_pres)) %>%
  mutate(sample="Full (n=1477)")

d2b <- datb[!is.na(datb$performance_health_pres),] %>% 
  group_by(treatment,performance_health_pres) %>% 
  summarise(count=n()) %>% 
  mutate(perc=count/sum(count)) %>%
  mutate(meanc=mean(datb[!is.na(datb$performance_health_pres) & datb$treatment == 0,]$performance_health_pres)) %>%
  mutate(meant=mean(datb[!is.na(datb$performance_health_pres) & datb$treatment == 1,]$performance_health_pres)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

pres_power = ggplot(d2, aes(x = factor(performance_health_pres), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  theme_bw() + 
  geom_text(aes(x=2, y=30, label=paste0("control:", "mu =","=", round(meanc, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=2, y=28, label=paste0("treatment:", "mu =","=", round(meant, 2))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("Very\ndissatisfied","Dissatisfied","Neither","Satisfied", "Very\nsatisfied")) + 
  labs(#title = "Perception of President's Performance in Health Service Provision",
       x = "Satisfaction with President's Performance",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "black"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.1, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/figure_4.tiff"), pres_power, width = 8, height = 5, dpi = 700)
ggsave(file = here("plots/pres_performance_health_wrap.pdf"), pres_power, width = 8, height = 5)

####################################################################################
# President general performance plot
d2 <- dat[!is.na(dat$performance_pres),] %>% 
  group_by(treatment,performance_pres) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(dat[!is.na(dat$performance_pres) & dat$treatment == 0,]$performance_pres)) %>%
  mutate(meant=mean(dat[!is.na(dat$performance_pres) & dat$treatment == 1,]$performance_pres)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Full (n=1477)")

d2b <- datb[!is.na(datb$performance_pres),] %>% 
  group_by(treatment,performance_pres) %>% 
  summarise(count=n()) %>% 
  mutate(meanc=mean(datb[!is.na(datb$performance_pres) & datb$treatment == 0,]$performance_pres)) %>%
  mutate(meant=mean(datb[!is.na(datb$performance_pres) & datb$treatment == 1,]$performance_pres)) %>%
  mutate(perc=count/sum(count)) %>%
  mutate(sample="Restricted (n=547)")

d2 = rbind(d2,d2b)

pres_power = ggplot(d2, aes(x = factor(performance_pres), y = perc*100, fill = factor(treatment))) +
  geom_bar(stat="identity", position = "dodge") +
  theme_bw() + 
  geom_text(aes(x=2, y=30, label=paste0("control:", "mu =","=", round(meanc, 1))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  geom_text(aes(x=2, y=28, label=paste0("treatment:", "mu =","=", round(meant, 1))),  size = 3, parse = T, check_overlap = TRUE) + #, family="CM Roman"
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(d2$perc*100)+1), labels=function(x) paste0(x,"%")) +
  scale_x_discrete(labels= c("Very\ndissatisfied","Dissatisfied","Neither","Satisfied", "Very\nsatisfied")) + 
  labs(#title = "Perception of President's Performance in General",
       x = "Satisfaction with President's Performance",
       y = NULL) + 
  scale_fill_manual(values=c("#999999", "black"),
                    labels=c("Control", "Treatment")) +
  theme(legend.title=element_blank(),
        plot.title = element_text(size=14,hjust = 0.5),
        #text=element_text(family="CM Roman"),
        legend.position = c(0.1, 0.9),
        legend.background = element_blank()) +
  facet_wrap(~sample)

ggsave(file = here("plots/figure_5.tiff"), pres_power, width = 8, height = 5, dpi = 700)
ggsave(file = here("plots/pres_performance_wrap.pdf"), pres_power, width = 8, height = 5)